%% it computes firm policy
 
mu=0;dd=100;
get_firm_policy_analyticalcase;
b_bar_anal=b_bar;
%% grid
if calib==1
    b_vec_d = 0.0:db:b_bar*2;
    N       = length(b_vec_d);
    %% matrices of derivatives
    Tb_f = zeros(N,N);
    for b=1:N-1
        Tb_f(b,b)   = -1/db;
        Tb_f(b,b+1) =  1/db;
    end
    Tb_f=sparse(Tb_f);

    Tb_b = zeros(N,N);
    for b=1:N-1
        Tb_b(b+1,b+1)=   1/db;
        Tb_b(b+1,b)  =  -1/db;
    end
    Tb_b = sparse(Tb_b);
    Tbb  = zeros(N,N);
    Tbb(1,1)      =  -1/db^2;
    Tbb(1,2)      =   1/db^2;
    Tbb(N,N)      =  -1/db^2;
    Tbb(N,N-1)    =   1/db^2;
    for b=2:N-1
        Tbb(b,b-1)=  1/db^2;
        Tbb(b,b)  = -2/db^2;
        Tbb(b,b+1)=  1/db^2;
    end
    Tbb=sparse(Tbb);
end

drift_b    = -rho*b_vec_d';
diff_b     =  sigma^2/2*(b_vec_d.^2)';
u_prime_l  = u1_l+u2_l*(1+gamma)*(min(1,b_vec_d/b_bar)).^gamma;
u_second_l = (gamma./b_vec_d).*(u_prime_l-u1_l);
u_third_l  = (gamma./b_vec_d).*u_second_l;
ell        = (-u_prime_l./u_second_l)*dr;ell(1)=max(ell(1),ell(2));
inde_b_bar = find(b_vec_d>b_bar,1)-1;
ell(b_vec_d>b_bar)=ell(inde_b_bar);
inde_db    = (ell./b_vec_d-rho)';inde_db(1)=100;

Tb                = Tb_f;
Tb(inde_db<=0,:)  = Tb_b(inde_db<=0,:);
Tb_db             = sparse(Tb.*(drift_b*ones(1,N)));
Tbb_dbb           = sparse(Tbb.*(diff_b*ones(1,N)));
AA                = Tb_db + Tbb_dbb;
%Check transition matrix conditions
if max(abs(sum(AA,2)))>10^(-10)
    disp('Improper Transition Matrix')
    return
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Analytical formulas
pp_v       = spline(b_vec,v_l);
v_anal     = -varphi*b_vec_d;
v_anal(b_vec_d<=b_bar)     = (ppval(pp_v,b_vec_d(b_vec_d<=b_bar)))';
x_anal     = -u_prime_l;
ell_anal   = ell;
inve       = zeros(N,1);
drift_I    = zeros(N,1);
inde_bi0   = 1:inde_b_bar;

%% model with investment
x_approx=1/100000;v_approx=1/100000;
a=.5; v_tresh = -(1-phi)*varphi*b_vec_d'+phi*ppval(pp_v,a*b_vec_d');v=ppval(pp_v,b_vec_d');

inde_db           = ones(N,1);
Tb                = Tb_f;
Tb(inde_db<=0,:)  = Tb_b(inde_db<=0,:);
options = optimset('Display','off');
tic
crit_v=1;inde_sim=0;
while crit_v>v_approx && inde_sim<50
    inde_sim=inde_sim+1;

    %obtain the investment given the value function
    drift_I_old      = drift_I;
    iota             = zeta*max(0,A*(v-(Tb*v).*b_vec_d')).^(1/theta);
    drift_I          = iota; drift_I(end) = 0;
    drift_b          =-(drift_I+rho).*b_vec_d';
    if zeta>0
        inve         = zeta^(-theta)/(1+theta)/A*(iota).^(1+theta); inve(end) = 0;
    else
        inve=zeros(N,1);
    end
    %solve for the new value function
    uu                = 1/nu-chi/A-(coupon+rho+varphi*delta)*b_vec_d' - inve;
    Tb_db             = sparse(Tb.*(drift_b*ones(1,N)));
    AA                = Tb_db + Tbb_dbb;
    get_LCP;v_new=v_lcp;inde=find(s_new<0.0001,1)-1;b_bar=b_vec_d(inde);
    if isempty(b_bar), b_bar=b_vec_d(end-1);end
    pp_v           = spline(b_vec_d,v_new); v_tresh = -(1-phi)*varphi*b_vec_d'+ phi*ppval(pp_v,a*b_vec_d');
    inde_b_bar     = find(b_vec_d>b_bar,1)-1;
    inde_bb        = inde_b_bar+1:N;
    inde_bi        = 1:inde_b_bar;

    %evaluate the guess
    crit_v     = max(abs(v_new(inde_bi)-v(inde_bi)))./mean(abs(v(inde_bi)));
    v          = v_new;

    %update the recovery rate
    u_prime_l  = (Tb*v);u_second_l = (Tbb*v);
    if inde_sim<500
        x_0         = -u_prime_l;x_0(inde_bb)=a_targ;    pp_x=spline(b_vec_d',x_0);
        [a ,fv]= fsolve(@(a)  a_targ - a*phi*ppval(pp_x,a*b_bar),a_targ/phi+0.1 ,options);
        a=min(0.7,a);
    end
end

ell        = (-u_prime_l./u_second_l)*dr;ell(1)=ell(2)*ell(2)/ell(3);ell(end)=ell(end-1)*ell(end-1)/ell(end-2);
x_0         = -u_prime_l;    pp_x=spline(b_vec_d',x_0);
x_0(inde_bb)=(1-phi)*varphi+phi*a*ppval(pp_x,a*b_vec_d(inde_bb));
x=x_0;

v_uncon0=v;
ell_uncon0=ell;
drift_I_uncon0=drift_I;
inde_bi_uncon=inde_bi;
x_uncon0=x;
inve_uncon0=inve;
dividend_uncon0         = 1/nu  +x_uncon0.*ell_uncon0 - (coupon+rho)*b_vec_d' - inve_uncon0;

%% model with  investment and adjustment frictions
% initial guess
inde_b_star= find(1/nu + x.*ell-(coupon+rho)*b_vec_d' - inve<0,1);
if isempty(inde_b_star)
    inde_b_star=N+1;
end

inde_db         = ones(N,1);
Tb              = Tb_f;
Tb(inde_db<=0,:)= Tb_b(inde_db<=0,:);

crit_x_out = 1;
tic
while crit_x_out>x_approx

    crit_v=1;inde_sim=0;
    while crit_v>v_approx && inde_sim<200
        inde_sim=inde_sim+1;
        v_prime         = Tb*v;pp_vprime=spline(b_vec_d(inde_bi)',v_prime(inde_bi));
        v_prime(inde_bb)=-(1-phi)*varphi+phi*a*ppval(pp_vprime,a*b_vec_d(inde_bb)');
        v_secon         = Tbb*v;v_secon(1)=v_secon(2)^2/v_secon(3);
        v_secon(inde_bb)= 0;

        %% compute x from HJB
        %drift  in b when constrained
        drift_db                 = ell+(sigma^2-rho-drift_I).*b_vec_d';
        %transition matrix
        Tb_db                    = sparse(Tb.*(drift_db*ones(1,N)));
        AA                       = Tb_db + Tbb_dbb;
        %cash flows
        uu                       = (coupon+rho+varphi*delta)*ones(N,1);
        %transition matrix when constrained:
        BB                       = (1/dd+r-dr+delta+rho)*speye(N)-AA;
        %transition matrix when unconstrained to the left of b*:
        BB(1:inde_b_star-1,:)    = zeros(inde_b_star-1,N);
        BB(1:inde_b_star-1,1:inde_b_star-1)    = speye(inde_b_star-1);
        %transition matrix after failure:
        BB(inde_bb,:)              = zeros(N-inde_b_bar,N);
        BB(inde_bb,inde_bb)        = speye(N-inde_b_bar);

        crit_x=1;x_new=x_0;x=x_0;dd=100;x_old=x;
        while crit_x>x_approx
            %cash flows
            bb                      = uu + x/dd;
            %terminal condition at point to the left ot bstar
            bb(1:inde_b_star-1)     = -v_prime(1:inde_b_star-1);
            %terminal condition at points to the right of bbar
            bb(inde_bb)             = -ppval(pp_vprime,a*b_vec_d(inde_bb)')*phi*a;
            %evaluate
            x_new              = BB\bb;
            x                  = x_new;
            crit_x             = max(abs(x_new-x));
            x = x_new;
        end

        %obtain implied lamdda in no-bankruptcy region
        lambda                    = zeros(N,1);
        lambda(inde_b_star:end)   = max(0,(-v_prime(inde_b_star:end))./x(inde_b_star:end)-1);

        %obtain investment policy
        drift_I_old      = drift_I;
        iota             = zeta*max(0,A*(v-(Tb*v).*b_vec_d')./(1+lambda)).^(1/theta);
        drift_I          = iota ;
        drift_I(inde_bb) = 0;
        inve             = zeta^(-theta)/(1+theta)/A*(iota).^(1+theta);
        inve(inde_bb)    = 0;

        %update ell
        ell_old = ell;
        ell(inde_b_star:end)=(inve(inde_b_star:end)+(coupon+rho)*b_vec_d(inde_b_star:end)'-1/nu)./x(inde_b_star:end);

        %solve for the value function
        %flow net return in unconstrained region
        uu                       = 1/nu-chi/A -(coupon+rho+varphi*delta)*b_vec_d' - inve;
        %flow net return in constrained region
        uu(inde_b_star:end)      = -chi/A-varphi*delta*b_vec_d(inde_b_star:end)';
        %drift in b in constrained region
        drift_b                  = (inve+(coupon+rho)*b_vec_d'-1/nu)./x -(rho+drift_I).*b_vec_d';
        %drift in b in unconstrained region
        drift_b(1:inde_b_star-1) =-(rho+drift_I(1:inde_b_star-1)).*b_vec_d(1:inde_b_star-1)';
        Tb_db             = sparse(Tb.*(drift_b*ones(1,N)));
        AA                = Tb_db + Tbb_dbb;
        %update value function and exit tresholds
        get_LCP;
        v_new= v_lcp;inde=find(s_new<0.000001,1)-1;b_bar =b_vec_d(inde);
        if isempty(b_bar), b_bar=b_vec_d(end-1);end
        pp_v           = spline(b_vec_d,v_new); v_tresh   = -(1-phi)*varphi*b_vec_d'+phi*ppval(pp_v,a*b_vec_d');
        inde_b_bar     = find(b_vec_d>b_bar,1)-1; inde_bb = inde_b_bar+1:N;%bankruptcy region
        inde_bi        = 1:inde_b_bar; %active region

        %check the inner loop guess
        crit_v     = max(abs(v_new-v))./mean(abs(v))+ max(abs(drift_I-drift_I_old));%+ max(abs(ell-ell_old))
        v          = v_new;
        %drift_I    = 0.5*drift_I+0.5*drift_I_old;
    end

    %% update ell and constrained region
    ell_uncon        = -v_prime./v_secon*dr;ell_uncon(inde_bb)=0;
    iota_uncon       = zeta*max(0,A*(v-(Tb*v).*b_vec_d')).^(1/theta);
    drift_I_uncon    = iota_uncon;
    inve_uncon       = zeta^(-theta)/(1+theta)/A*(drift_I_uncon).^(1+theta);
    dividend         = 1/nu  -v_prime.*ell_uncon - (coupon+rho)*b_vec_d' - inve_uncon;
    inde_b_star_old  = inde_b_star;
    inde_b_star_cand = find(dividend<0,1);
    if ~isempty(inde_b_star_cand)
        if inde_sim>5%update only if distance from zero is larger than at first negative
            if dividend(inde_b_star)>abs(dividend(inde_b_star_cand))
                inde_b_star      = inde_b_star_cand;
            end
        else
            inde_b_star=inde_b_star_cand;
        end
    else
        inde_b_star=N+1;
    end

    %evaluate outer loop guess
    crit_x_out = abs(inde_b_star-inde_b_star_old)
    %update ell at new b_star
    ell(1:inde_b_star-1)        = ell_uncon(1:inde_b_star-1);
    ell(inde_b_star:end)        = 1./(x(inde_b_star:end)).*((coupon+rho)*b_vec_d(inde_b_star:end)' + inve(inde_b_star:end)- 1/nu);
end
toc


if crit_v>v_approx,warning('v not converged'),end
if b_bar==b_vec_d(end-1);warning('never fails'),end

pp_ell = spline((b_vec_d(inde_bi(1:end))),ell(inde_bi(1:end)));
pp_inv = spline((b_vec_d(inde_bi(1:end))),inve(inde_bi(1:end)));
pp_drift = spline((b_vec_d(inde_bi(1:end))),drift_I(inde_bi(1:end)));
pp_x   = spline((b_vec_d(inde_bi(1:end))),x(inde_bi(1:end)));
pp_v   = spline((b_vec_d(inde_bi(1:end))),v(inde_bi(1:end)));

b_star = b_vec_d(inde_b_star);